clear; close all; clc;

%% 1. 模型参数设置
alpha   = 4.5;        % (α, σ, μ) = (4.5, 0.14, 0.001)
sigma   = 0.14;
mu      = 0.001;
k       = 0.1;        % 耦合强度
phi0    = -1;         % 忆阻器初始状态
epsilon = 1;          % phi 的更新系数

idx = 0:3500;

%% 2. 初始状态矩阵
% 每行：[X1_0, φ0]，其中 X2_0 取 -X1_0，Y1 和 Y2 均固定为 0
IS = [ 0.1, -1;
      -0.1,  1;
       0.1,  1];

%% 2. 迭代设置
N_total   = 3500;      % 总迭代步数
% 对于绘图而言，我们希望获取完整的时间序列，因此在调用 Rulkov.m 时设置 N_cut = 0
N_cut     = 2500;      % 不丢弃瞬态

%% 3. 调用 Rulkov.m 函数获得激发序列
% X1_save 和 X2_save 分别为两个神经元的激发变量序列
for i = 1:size(IS,1)
    k      = IS(i,1);
    phi0  = IS(i,2);
    
    % 调用修改后的 Rulkov 函数（支持自定义初始 X1_0）
    [X1_save, X2_save, ~] = Rulkov(alpha, sigma, mu, k, phi0, N_total, epsilon, N_cut);
    
    % 绘图
    subplot(3,1,i);
    idx = N_cut+1 : N_total;
    plot(idx, X1_save(idx), 'r-' , 'LineWidth',1); hold on;
    plot(idx, X2_save(idx), 'b--', 'LineWidth',1);
    xlabel('n');
    ylabel('X_1, X_2');
    title(sprintf('初始: k=%.2f, \\phi_0=%.2f', k, phi0));
    legend('X_1','X_2','Location','Best');
end